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ABSTRACT 

We present axisymmetric simulations of the coupled convective and radiative regions in the Sun in 
order to investigate the angular momentum evolution of the radiative interior. Both hydrodynamic 
and magnetohydrodynamic models were run. We find an initial rapid adjustment in which the dif- 
ferential rotation of the convection zone viscously spreads into the radiative interior, thus forming a 
"tachocline" . In polar regions the subsequent spread of the tachocline is halted by a counter-rotating 
meridional circulation cell which develops in the tachocline. Near the equator such a counter-rotating 
cell is more intermittent and the tachocline penetration depth continues to increase, albeit more slowly 
than previously predicted. In the magnetic models we impose a dipolar field initially confined to the 
radiative interior. The behavior of the magnetic models is very similar to their non-magnetic counter- 
parts. Despite being connected to the convection zone, very little angular momentum is transferred 
between the convective and radiative regions. Therefore, while it appears that a magnetic field is not 
necessary to stop the tachocline spread, it also does not promote such a spread if connected to the 
convection zone. 

Subject headings: solar physics, magnetohydrodynamics 



1. INTRODUCTION 



Helioseism ic observations ([Brown eFaH [1981 



iGoode et al.l I1991D have revealed a surprising pic- 
ture of the solar rotation profile. These observations 
show that the differential rotation observed at the 
surface, with the equator spinning faster than the poles, 
persists to the base of the convection zone with very 
little radial variation. On the other hand the radiative 
interior is revealed to be rotating uniformly. The 
transition between these very different rotation profiles 
is now commonly known as th e tachocline, with a stil l 
undetermined extent, although lElliott fc Goughl ([1999f ) 
estimate a width of O.O19i?0. 

The observed rotation profile is puzzling for many rea- 
sons. First, it was previously thought that the rota- 
tion profile of the solar convection zone would obey the 
Taylor-Proudman (TP) constraint, in which rotation is 
constant on cylinders. Indeed, it has proven quite diffi- 
cult to break this constraint and recover differential ro- 
tation constant on cones, similar to the observed pro- 
file. [Rempcl (2005) showed that the TP constraint could 
be overcome with a latitudinal entropy gradient imposed 
in a subadiabatic tachocline. This thermal wind would 
break the TP constraint and lead to more realistic dif- 
ferential rotation in the con vection zone. F ollow up nu- 
merical experiments by Mie sch et al.l ([2006[ ) showed that 
indeed, when a latitudinal temperature gradient of order 
lOK is imposed as a boundary condition at the base of 
their model convection zone, solar-like differential rota- 
tion is recovered. 

The studies of the convection zone differential rota- 
tion discussed above suggest the importance of the un- 
derlying tachocline and adjacent stable region in un- 
derstanding the differential rotation of the convection 
zone. This leads to the second (and third) puzzling 
aspect of the solar rotation profile. Despite the con- 
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stant forcing by the overlying differential rotation of the 
convection zone the tachocline remains thin and the ra- 
diative interior is maintained in uniform rotation. Al- 
though a thin tachocline does not guarantee a uniformly 
rotating radiative interior, the two questions have of- 
ten been t reated as one. I n the defining work on the 
tachocline, iSpiegel fc ZahnI ([1992t ) attribute the thinness 
of the tachocline to anisotropic turbulence. They ar- 
gue that, similar to the Earth's atmosphere, the solar 
tachocline, lying mostly in the stably stratified interior, 
opposes vertical turbulence. They show that the width 
of the tachocline is then determined by the degree of 
anisotropy in the turbulence. Simulations carried out by 
lElliottl ([1997f ) confirmed this trend. 

An alternative sol ution was offered by 
IGough fc Mclntvrd ((l998h . They argue that there 
are many examples in the Earth's atmosphere in which 
two-dimensiona. 1 turbu lence of the sort envisioned by 
ISpiegel fc ZahnI ([1991 conserve potential vorticity and 
thus drive a system away from uniform rotation. These 
authors take the helioseismic results as an indication 
that there must exist a large scale magnetic field in the 
radiative interior. Such a field would easily establish 
uniform rotation of the r adiat ive interior via Ferraro's 
isorotation law ([Ferrarol Il937| ) and a thin tachocline 
could be formed as something of a magnetic boundary 
layer. 

The magnetic solution to the rotation pro- 
file of the radia tive interior was not unique to 
IGough fc Mclntvri ([T998 ). it was investigated 
by other authors (jC harbo nneau fc MacGregori 

19931 iRudiger fc K itchatinol 119971 : 

MacGregor fc CharbonneaiJ 1999). In particular 
MacGregor fc Charbonneaul (|1999l) elucidated a problem 



with this theory that is now known as the "confinement" 
problem. That is, magnetic field naturally diffuses 
and any remnant field in the radiative interior would 
diffuse into the convection zone. Once the field is 
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in contact with the differentiahy rotating convection 
zone, a toroidal field and associated Lorentz force act 
to spread the differential rotation of the convection 
zone into the radi at ive in terior. The contribution of 
iGough fc Mclntvr^ (|1998f l was to suggest that this 
diffusion could be controlled by meridional circulation 
driven in the convection zone and burrowing into the 
tachcoline; effectively holding down the magnetic field. 
Subsequent studies on the ability of a meridional fiow 
to confine the tachocline h ave produced mixe d results: 
iGaraud (2002); Kitchatin ov fc Riidigeil (|2006[ ) initially 
showed that such confinement by meridional circulation 
was possible. Time-depende n t, thr ee-dimensional (3D) 
simulations bv iBrun fc ZahnI (|2006[ ) indicated that that 
it was not. The differences in these works amounted to 
how deep the meridional circulation, driven in the con- 
vection zone, could penet rate into the radiative interior , 
and with what strength. iKitchatinov fc Riidigerl ()2006l ) 
estimated that a flow speed of 10^ cm/s penetrating 
lO^cm would be sufficient to confine the field. 

The strength of meridional circulation and it's ability 
to penetrate the stable tachocline depends, in turn, on 
how the boundary condition with the convection zone is 
treated. In all of these studies only the radiative interior 
was modeled, with the convection zone treated purely 
as a boun dary condition which impos ed the differential 
rotation. IGaraud fc Brummelll ()2008f ) showed that the 
depth of penetration of meridional circulation into the 
radiative interior and its strength depended sensitively 
on the boundary conditions assumed. In reality, the 
convection zone imposes a complex boundary condition 
with (at least) large scale infiows and outfiows, penetra- 
tion of plumes and pumping of magnetic field. However, 
coupling the convection zone and radiative interior is a 
difficult undertaking because of the vast discrepancy in 
the governing timescales for the convective and radiative 
zones in the Sun. Here we report on axisymmetric, time- 
dependent simulations which solve self-consistently both 
the convection zone and the underlying radiative region. 

2. NUMERICAL MODEL 

2.1. Equations 

The numerical method we em ploy here is identical to 
that described in IRoger^ (|2011[ ). here we give a brief re- 
view and refer the reader there for more details. We 
solve the full set of axisymmetric MHD equations in the 
anelastic approximation: 



— = V X (u X B) + r/V^B 

ot 



(5) 



V • (pu) = 0. 
V -B = 



(1) 
(2) 



(9ii 

— + (u • V)u = -VP - Cgf + 2(u X n) 

-fi(JxB) + i7(v2u + iv(V-u)) (3) 
P 3 



dT 

- + (v.V)r^ 



dT — 

- (7 - ^)Thp) 



Equations 1 and 2 ensure the mass and magnetic flux 
are conserved. In Equation 3, the momentum equation , 
C denotes the co-density ([Braginskv fc RobertsI 119951 : 
[Rogers fc Glatzmaieri I2005[) . g the gravity, fl the rota- 
tion rate, v and k the viscous and thermal diffusivity 
respectively. The inverse density and thermal diffusivity 
scale heights are denoted by hp and /i^, 7 is the adiabatic 
index and all other variables take their usual meanings. 
In all of the above equations reference state variables 
are denoted by an overbar. Equation 4 reresents the en- 
ergy equation written as a temperature equation, where 
Ur represents the radial velocity. The first term on the 
right hand side (rhs) of equation 4 represents the differ- 
ence between the reference state temperature gradient 
and the adiabatic temperature gradient and allows us 
to represent strongly subadiabatic regions in addition to 
convection zones. Equation 5 is the magnetic induction 
equation, in which rj represents the magnetic diffusivity. 
The diffusivities have the following forms: 



P 

const 



(6) 



+ (7 - l)ThpUr + jk[V'T + {hp + K)^] (4) 



where Km is 10 , is 5 x 10 and 77 is 5 x 10 , giving 
units of cm} js. Values at the base of the convection 
zone are shown in Table 1. It is worth mentioning here 
that because our thermal diffusivity, k is 10^ larger than 
solar the flux through the simulation is 10^ x solar. Using 
Mixing Length Theory (MLT) this would imply velocities 
which are (10^)^/3 ^ ^qOx solar. 

These equations are solved in axisymmetric spherical 
coordinates using a spherical harmonic decomposition in 
latitude and a finite difference scheme in the radial di- 
rection. The maximum spherical harmonic degree, Imax, 
is 340, so that the grid resolution is 512 zones in lati- 
tude. The radial resolution is 1500 zones, with 600 zones 
dedicated to the radiative interior and 900 zones cov- 
ering the convection zone. We solve the above equa- 
tions from 0.15 R© to 0.95 R©. The radiative region 
extends from 0.15 R© to 0.71 R© and the convection 
zone extends from 0.71 R© to 0.90 R©, for numerical 
reasons we impose an additional stably stratified region 
from 0.90 R© to 0.95 R©. The reference state thermo- 
dynamic variables are taken from a polynomi al fit to the 
standard solar model from Christcnscn-Dals gaard et al.l 
(1996). We should note that because of the convection 
and the gravity waves in the deep radiative interior our 
timestep is severely limited and therefore, we are unable 
to run these simulations very long. Our average timestep 
is approximately 10s and all of these simulations have run 
« 2 X lO^s. 

2.2. Model Setup 

We recognize that without fully 3D simulations we 
would be unable to recover the observed solar differential 
rotation profile in the convection zone. However, we are 
mainly interested in studying the response of the under- 
lying radiative interior to this differential rotation in the 
presence of physical phenomena such as convecti ve over- 
shoot and magnetic pumping. Therefore, as in iRogerj 
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Symbol 


Sun 


Simulation 


thermal diffusivity 


Kj 


10^ 


1.1 X 10^^ 


magnetic diffusivity 


Ti 


10^ 


5 X 10^^ 


viscous diffusivity 


y 


30 


2.5 X 10^^ 


Praudtl number 




3 X lO^*' 


0.22 


magnetic Prandtl nb 


u /r) 


3 X IQ-'^ 


0.5 


Froude number 


{n/Nf 


2 X 10-6 


2 X 10-s 


Ekman number 




2 X 10-15 


2 X lO-'' 


rotation frequency 


n 


2.7 X lO-'^ 


2.7 X lO-'* 


Ohmic diffusion time 


tB 


2.5 X 10^* 


5.0 X 10** 


viscous diffusion time 




8.333 X 10" 


109 


Eddington-Sweet time 




2.8 X lO^s 


2.8 X 1013 



TABLE 1 

Values for various parameters in the simulations. The 

VALUES listed FOR THE DIFFUSIVITIES ARE THOSE AT THE TOP OF 
THE RADIATION ZONE AND HAVE UNITS cni^/s. LIKEWISE, THE 

timescales quoted are simply calculated using the 
diffusivities at the top of the radiative interior. 
Diffusion times are calculated using the depth of the 
radiation zone and are quoted in seconds. each model was 
run approximately 2 x 10® s and required approximately a 
few hundred thousand processor hours. 



(|2011f ) we artificially force the differential rotation of the 
convection zone to be that inferred from helioseismology 
((Thompson et al. 1996) and described by: 



nc,^A + B cos^ e + C cos"* I 



(7) 



where A is 45 6nHz, B is -72nHz and C is -42nHz. Unlike 
IRoeersI ((20T1 we force the differential rotation only in 
the convection zone, with a step function to the uniform 
rotation of the radiative interior^ The differential rota- 
tion is imposed through a forcing term on the azimuthal 
component of the momentum equation, which is present 
only in the convection zone: 



d_ 

dt \r sm ( 



U4, 



(V • uu^)^ [(V X B) X B]^ (2u X 



r sm 



piJLr sm ( 
vu^ _ 

r-3 sin"^ 9 



F 



(8) 



where the last term, F, represents the forcing function, 
defined as: 



F 
F 



si'(r, e) 



r sin 6 



0. 



for r > O.71i?0 
for r < 0.71i?o 



(9) 



where ^l'{r,9) represents the observed differential ro 
tation relative to the constant n and r is a forcin; 
timescale. In all simulations ri=441nHz and r is 10'*s 

We ran three models. One purely hydrodynamic model 
and two magnetic models. In the magnetic models we 
start with a hydrodynamic model until convection is suf- 
ficiently established at which point we imposed a dipolar 
field with the form: 

1 While we have run models both with an imposed tachocline 
and those without, we feel that imposing a tachocline does not 
accurately represent the situation we are trying to investigate. This 
is because if a tachocline is forced, there is no stress presented to 
the radiative interior. 

^ We have run models with t = 10^ and 10'^ for limited time 
(not as long as the full simulations presented here, but much longer 
than the forcing time) and see very little difference. 



0.71R 







(10) 



so that initially all field lines close just beneath the con- 
vection zone. Figure [T] shows a sch ematic o f the initial 
setup and model, reproduced from iRogersI (|2011[ ). Al- 
though we ran models with two magnetic field strengths. 
Bo = 40G and Bo = 4000G, all of the results presented 
will be only for the strong field case because the weak 
field case is virtually identical to the non-magnetic case. 

In all models the boundary conditions are stress-free 
and impermeable at the top and bottom boundaries. The 
temperature boundary conditions are constant tempera- 
ture perturbation on top and bottom and the magnetic 
field is matched to an internal potential field at the bot- 
tom of the domain and an external potential field at the 
top of the domain. The current density, J/ is zero on the 
top and bottom boundary. 

3. RESULTS 

3.1. Meridional Circulation 

The fundame n tal assumption of the 

iGough &: Mclntvri ()1998f ) model is that meridional 
circulation driven by the differential rotation in the 
convection zone is able to penetrate (at least minimally) 
the radiative interior and confine the magnetic field. 
One of the persistent issues in testing this model has 
been regarding how deep this meridion al circulation pen- 
etrates the radiative interior, if at all. iGilman fc MiescH 
(2004!) argued that the depth of penetration of the 
meridional circulation driven within the convection 
zone was limited to an Ekman depth and therefore, 
would penetrate no further than the overshoot layer. 
Garaud & Brummell (2008) showed that the circulation 
driven beneath the convection zone depends strongly 
on the boundary conditions imposed at the top of the 
radiative layer and that radial forcing at the interface 
could produce small, but sustained fiows deep within 
the interior. 

The simulations presented here couple the convective 
and radiative regions. Moreover, the radiative interior 
is sufficiently "stiff" , using the Brunt- Vaisala frequency 
from the standard solar model. Therefore, these simula- 
tions are able to address the penetration of meridional 
circulation in the most self-consistent manner to date. 
Furthermore, because we artificially impose the differen- 
tial rotation in the convection zone we expect a faithful 
representation of the meridional circulation driven in the 
convection zone from the differential rotation. 

In Figure [2] we show the latitudinal velocity, both in 
time-average and time snapshots. In this figrue red rep- 
resents flow towards the south pole, while blue repre- 
sents flow towards the north pole. One can immediately 
see the poleward flow at the surface, with equatorward 
flow near the base of the convection zone (the overlaid 
red line in (a) and (b) shows the base of the convection 
zone). In both the magnetic and non-magnetic cases the 
time-averaged meridional circulation shows more struc- 
ture near the equator and is not well described as a single 
cell, but this is likely due to needing longer time averages. 
The amplitudes of the meridional circulation within the 
convection zone show velocities of approximately 9 x 10^ 
cm/s at the top of the convection zone and 5 x 10'' cm/s 
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near the base of the convection zoneO One can also see 
hints of internal gravity waves (IGW), particularly in the 
non-magnetic case and in the time-snapshots (these can 
be seen more clearly in Figure |3|). 

While Figure [5] clearly shows the general circulation 
set up in the convection zone it is hard to discern the 
amplitudes of the meridional flows within the radiative 
interior. In Figure [3] we show these flows more clearly. 
In that figure we show the absolute magnitude of the 
both the radial and latitudinal velocities as a function of 
radius below the convection zone for three different lat- 
itudes. Very little latitudinal dependence is seen. How- 
ever, one can clearly see a rapid decrease in depth of the 
flow velocities, which are down nearly four to five orders 
of magnitude by O.GOi?©. Figure [3^ and [3]: show this 
amplitude for non-magnetic cases while [Sb and[3ji show 
it for the magnetic cases. One can see a slightly faster 
decay in the magnetic cases as well as shorter vertical 
wavelengths of the mixed gravity-alfven waves, which is 
expect ed given the dispersion rela tion for these mixed 
waves (jMacGregor fc Rogersl[20Tl . 

Accordi ng to the analyt ic expression derived by 
ICaraud fc Acevedo-Arreguinl ([2009 ). the amplitude of 
the radial velocity should fall off exponentially with a 
length-scale determined by the parameter a = y/PrN/fl, 
where N is the Brunt- Vaisala frequency. In these simula- 
tions the Prandtl number and Brunt- Vaisala frequency, 
N, vary substantially over the domain, so a direct com- 
parison is difficult, but using N « lOOfiHz and a 
Prandtl number of 0.05 (the values just beneath the con- 
vection zone), the analytically predicted length-scale is 
« 9 X lO^cm. Investigation of the data plotted in Fig- 
ure [3] gives a length-scale of w 4 x lO^cm, in fair agree- 
ment with the analytic prediction (given the simplified 
comparison). 

By forcing the differential rotation in the convection 
zone we are artificially breaking t he Taylo r-Proudman 
constrain t there. As p redicted by iRempell ([2005) and 
tested by iMiesch et al.l (2006) in three-dimensional sim- 
ulations, breaking the Taylor-Proudman constraint could 
be achieved by a latitudinal temperature gradient within 
the tachocline (or as a bottom boundary condition in the 
case of the simulations). Here we can ask the question 
in reverse; given the imposed differential rotation, is a 
latitudinal temperature gradient recovered? In Figure U 
we show the temperature fiuctuation as a function of lat- 
itude at three radii within the tachocline (the tachocline 
will be defined in the next section, see figure caption for 
actual depths beneath the convection zone) . Here we see 
poles which are much warmer than the equator at the 
top of the tachocline (solid line) and slightly lower in 
the tachocline (dotted line). However, we also see a fair 
amount of variation at lower latitudes. In addition, near 
the middle of the tachocline (dashed line) this temper- 
ature variation in latitude switches sign (particularly at 
high latitudes) , so that the poles are slightly cooler than 
the equator. 

^ The amplitudes of the meridional circulation here are larger 
than those observed in the solar convection zone because of an 
increased thermal diffusivity, as discussed in Section 2.1. While 
this is not ideal, our overestimate of the flow speeds likely leads us 
to enhanced penetration of flows. Weaker flows would lead to less 
penetration. 



Investigation of Figure [2] indicates that the middle 
tachocline (where the latitudinal temperature gradient 
changes sign) is also assocated with a meridional circu- 
lation cell which rotates counter to the main meridional 
circulation driven in the convection zone. This counter 
rotating cell (CRC) is robust in the polar regions. Near 
the equator the cell is more time dependent, but longer 
time averages show the CRC becoming more prevalent. 
We should n ote that such a CRC has been seen in 3D sim- 
ulations by IMiesch et al.l (|2010f ) , however in that work 
they attribute this cell to their artificial impermeable 
boundary at the base of their convection zone. We con- 
firm here the presence of such a cell even under realistic 
boundary conditions. This CRC has important implica- 
tions for the spread of the tachocline as we will discuss 
in Section 3.4. 

3.2. Tachocline Penetration Depth 

In the simulations discussed here we do not impose a 
tachocline, but instead want to investigate the spread of 
differential rotation from the convection zone into the 
radiative interior. Thus, look at the formation of the 
tachocline and its evolution in time. Figure [S] shows both 
time-averaged and time snapshots of the azimuthal ve- 
locity. One can clearly see that the differential rotation 
quickly penetrates into the radiative interior. However, 
between Figure [Sis and [Sf, there is very little change in 
that penetration, although panel [5f represents a substan- 
tial increase in time. One way to measure the tachocline 
depth is by looking at the region of strongest radial 
shear in the azimuthal velocity. Looking at \d{v^/r)/dr\ 
as a function of radius and time we can then measure 
the tahcocline penetration depth as the depth beneath 
the convection zone at which this is half its peak value. 
This is shown in Figure [SI which shows the radial gra- 
dient of angular velocity as a function of radius at both 
high (left-hand panels) and low (right-hand panels) lat- 
itudes at various times. The dashed lines represent 
the secondary peak and the solid lines represent half of 
that peak. If the tachocline penetration depth is mea- 
sured as the distance between the base of the convection 
zone and the solid line, the average tachocline depth is 
O.O21i?0 (1.46 X lO^cm) at high latitudes and 0.022i?o 
(1.53 X lO^cm) at low latitudes0 Measured this way, the 
tachocline penetration depths we measure are consistent 
with helioseismic results. 

Another way to define tachocline thickness at a given 
latitude is the depth beneath the convection zone at 
which the angular velocity perturbation (initally zero) is 
X% of its forced value at that latitude within the convec- 
tion zone. We plot this value for X=50% in Figure[71 and 
note that the general behavior in time at different values 
of X is the same even if the amplitude is not. Figure [7| 
shows the tachocline thickness, measured this way, as 
a function of time, averaged over polar regions (0°-30°) 
and equatorial regions (60°-90°) for the non-magnetic 
case (sold hue) and the strong field case (dashed line). 
Also shown in that figure is the time derivative of the 
tachocline penetration depth as a function of time (pan- 
els (b) and (d), and zoomed in over the last half of the 

This provides an upper limit for this way to measure the 
tachocline depth. If we had measured from the secondary peak 
to its half max, the depth would be much smaller. 
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simulation in the insets). The first obvious feature of 
these results is that the tachocline thickness increases 
rapidly, and subsequent evolution is slow. In fact, in 
the polar regions the tachocline penetration depth has 
stopped increasing, as can be seen in Figure [7^, as well as 
[7t), where the time derivative is sometimes negative. In 
the equatorial regions the tachocline penetration depth 
is still increasing, however at a very slow rate, approx- 
imately 1-5 cm/s (as seen in the insets of panels b and 
d of Figure [7]), which would take longer than the so- 
lar age to spread through the entire solar radiative inte- 
rior. Measuring the tachocline depth this way we recover 
tachocline depths of w 8x 10^-10i°cm or O.12-O.16i?0. 

The tachocline widths measured this way are larger 
and inconsistent with helioseismic results. However, 
given the diffusivitics used in these simulations a thicker 
tachocline is not surprising. Furthermore, the most im- 
portant feature is that the tachocline thickness has vir- 
tually stopped spreading (particularly at high latitudes). 
At high latitudes, the tachocline thickness at the end of 
the simulation is 97% that at half the simulation time for 
the non-magnetic case and 98% for the magnetic case. At 
low latitudes the tachocline thickness is 2% higher at the 
end of the simulation compared to half of the simulation 
time for both magnetic and non-magnetic cases. 

The viscous spread of the tachocline should give a 
tachochne thickness A « R{t/U)''-/'^ . After 4 x lO'^s 
(the time of initial rapid tachocline spread in these sim- 
ulations) and using the viscous diffusivity at the top of 
the tachocline, the viscous depth is approximately 0.16 
i?0 (lO^^cm), using the viscous diffusivity at the bot- 
tom of the tachocline this depth is approximately O.OSRq 
(5.5 X lO^cm), in excellent agreement with the tachocline 
thicknesses seen her^l. However, if we were to assume 
viscous spread over the entire simulation time, we would 
get tachocline penetration approximately equal to the 
depth of the radiative interior. 

If, on the other hand, the subsequent evolution of 
the tachocline (after the rapid initial spread) obeyed 
the circulation-enhanced, "hyper-diffusive" spreading de- 
scribed in lSpiegel fc Zahiil (|l992) in which A cx t^^^, then 
the tachocline penetration depth would be approximately 
1.7 times larger at the end of the simulation than it is at 
1/3 of the simulation time. This is also not seen, clearly 
indicating that something is preventing both the viscous 
and the thermal hyper-diffusion spread of the tachocline, 
even in the non-magnetic casfl 

In the following sections, we refer to the "tachocline" 
in these models as defined by the radial gradient of the 
azimuthal velocity or O.O2i?0. Although we recognize 
that there are other ways to define the tachocline depth 
(as studied above) we find that the relevant force balance 
which stops the tachocline spread occurs in this narrow 
depth. We now turn to the question of what is stopping 
the spread of the tachocline. 

^ The init i al ra pid adiabatic response predicted by 
ISpiegel fc ZahnI 1)19921 ') occurs on a timescale of approximately 
2 X 10®s, shorter than what we see here. Furthermore, the depth 
over which it extends is predicted to be approximately 4 X lO^cm, 
again shorter than we see here. 

® In these simulations the Prandtl number w/k in the tachocline 
is approximately 0.05, but varies with depth. This is not low 
enough to ensure that the spread due to thermal hyper-diffusion is 
larger than that due to viscous diffusion, as in the Sun. 



3.3. The Spiegel & Zahn Model 

The defining work on the tachocline bv lSpiegel fc ZahnI 
(|1992') proposed that the spread of the tachocline could 
be stopped by an anisotropic turbulence, akin to that 
seen in the Earth's atmosphere. This is argued based 
on the fact that the bulk of the tachocline lies within 
the radiative interior, where the Brunt- Vaisala frequency 
is large and likely inhibits substantial radial motion. 
Since our non-magnetic model also shows a halt of the 
tachocline penetration depth this seems a likely candi- 
date. We can test this proposal by investigating the 
stresses within the tachocline. The Reynolds stress can 
be broken into mean and fluctuating components, such 
that: 

V • (UU^) = V • {UU'^) + V • (UM^) (11) 

where iT^ represents a time average and u'^ represents the 
fluctuation, so that = + u'^. 

The first term on the rhs of Equation 11, the divergence 
of the Reynolds stress, can be split into horizontal and 
radial terms and we can ask if there is any anisotropy. 



1 d{u'^ug sin 9) u'^ue cot ( 



(12) 

r sin 9 do r 

The anisotropy envisioned by ISpiegel fc ZahnI (|1992D 
would then require that the ratio of the radial to lati- 
tudinal terms on the rhs of Equation 12 be small. We 
show this ratio in Figure [8] as a function of latitude 
at various depths within the tachocline. In both mod- 
els and at all radii evaluated in the tachocline there 
is no detectabl e aniso tropy of the kind envisioned in 
ISpiegel fc ZahnI (| 19921) . This could be due to our lack 
of resolution, so we can not rule this model out com- 
pletely, but we can say it is not responsible for the halt 
of tachocline spread seen here. In fact, in much of the 
tachocline these Reynolds stresses are small compared to 
other terms in Equation 8 (see next section). 

3.4. Force Balance 

The above discussion then implies that the 
ISpiegel fc ZahnI ()1992[ ) model is not at work in these 
simulations. Furthermore, given that the tachocline 
spread is halted in the non-magnetic cases, we also 
see that a magnetic field is not necessary in stopping 
the tachocline spreacQ. This leads to the important 
question of what is responsible for limiting the tachocline 
penetration depth. For an answer we can look directly 
at the force balance within the tachocline to determine 
which forces are responsible for stopping this spread. 
The relevant equation is the azimuthal component of 
the momentum equation given by: 



d 



dt Vrsin6' 



1 ^ / N [(V X B) X 



+ 



r sm 

(2u X ri)^ 



p^,r sm f 



where F has been omitted because there is no forcing in 
this region and, as in Equation 11, the total azimuthal 

^ Although see magnetic section for further discussion. 
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velocity, is made up of the mean (time-averaged, u^) 
plus the fluctuating component, u'^. 

In Figure [9] we show the various force terms as a 
function of latitude, at two different radii within the 
tachocline, with the top four panels representing the non- 
magnetic case and the bottom four panels representing 
the strong magnetic cas^l- The left-hand panels repre- 
sent high latitudes, while the right hand panels represent 
low latitudes. In the following discussion we will con- 
centrate on the non-magnetic case, but note here that 
the magnetic case shows the same behavior. At high 
latitudes, at the top of the tachocline (top left panel) 
the total force (denoted by the solid line) is very small. 
The force there is dominated by the the Coriolis force 
(dashed line, defined as the sum of the first and fourth 
terms on the rhs of Equation 13) and the divergence of 
the Reynolds stress (dotted line, defined as the second 
term on the rhs of Equation 13). The divergence of the 
Reynolds stress in this region is largely due to convective 
overshoot and is therefore, unsurprisingly large. The sign 
of the Coriolis force at the top of the tachocline is nega- 
tive indicating its tendency to slow the poles. However, 
this Coriolis term is offset entirely by the Reynolds stress 
and diffusive terms, leading to very little net angular mo- 
mentum transport there. 

At the bottom of the tachocline, again at high lat- 
itudes, the Coriolis term, which is still dominant has 
changed sign and leads to an acceleration of the az- 
imuthal flow at high latitudes. The change of sign of the 
Coriolis force is because, at the bottom of the tachocline, 
the Coriolis force is acting on the counter-rotating merid- 
ional circulation cell discussed in Section 3.1, which is 
poleward. Therefore, this counter rotating cell acts to 
change the sign of the dominant force in the tachocline, 
effectively turning off the spread of the tachocline in ra- 
dius. 

At lower latitudes, the situation is more complicated. 
There is substantially more time and latitudinal varia- 
tion. This is because a persistent CRC has not developed 
at lower latitudes. We believe that over longer time aver- 
ages the low-latitude CRC will become more persistent 
and halt the spread there as well. This expectation is 
supported by Figure [7] which shows the spread of the 
tachocline nearing zero, even at low latitudes. 

3.5. Magnetic Models 

As stated above we ran three models: two with a mag- 
netic field and one without. Although we ran two mag- 
netic cases, as stated previously we only present results 
from our strong field case, for ease of presentation, but 
also because the results in the weak field case are vir- 
tually identical to the non-magnetic case. The results 
for the tachocline penetratio n depth shown in Figure [71 
and Figure m indicate that the lGough fc Mclntvrd (|l998f ) 
model can not be solely responsible for halting the spread 
of the tachocline. However, it is also expected that an 
internal poloidal field once "connected" to the differen- 
tially rotating convection zone would cause the differen- 
tial rotation to be spread into the radiative interior on 
approximately an Alfven crossing time. In Figure fTOl we 

* The magnetic stresses are hardly seen because they are sub- 
stantially lower than the other stresses in the tachocline, see section 
3.5. 



show the evolution of the poloidal (and toroidal) field, 
where we can clearly see that the poloidal field connects 
to the convection zone. Yet, despite this connection we 
see no additional spread of the tachocline in the magnetic 
case. 

For our strong field case, with an initial internal field 
strength of 4kG, assuming no diffusion, nor any turbulent 
diffusion due to small scale motion, the Alfven velocity 
varies from a few cm/s at the base of the convection 
zone to a few hundred cm/s in the deep radiative inte- 
rior. After the initial rapid spread of the tachocline as 
seen in Figure [71 that is at ~ O.6i?0, the Alfven velocity 
(va) is approximately 100 cm/s. Therefore, in the subse- 
quent lO^s after the initial spread of differential rotation, 
one would (very simply) expect an additional spread of 
Vat « lO^^cm if that spread were being mediated by a 
magnetic field. However, we see no spread of this mag- 
nitude. Indeed, we see very little effect of the magnetic 
field at all. 

There are (at least) three reasons for this (that this 
author can think of). The first, and most obvious, is 
that these simulations have not run long enough to see 
the additional effects of the magnetic field. The esti- 
mate above, is just that, a very crude estimate. Because 
we wanted our magnetic Prandtl number to be less than 
one, our magnetic diffusivity is rather large. Compound- 
ing this problem we have no dynamo and therefore, the 
field strength and Alfven velocity are decreasing (after 
an initial rise, see Figure [TT|) . Estimating the decay 
of the magnetic field strength using the imposed diffu- 
sivity, the amplitude of the magnetic field (and hence, 
Alfven velocity) would have dropped by approximately 
30% and hence, the spread might have only amounted to 
7 X lO^cm, but this would still be measurable. The pic- 
ture is more complicated though, because the evolution 
of the magnetic field is not governed solely by diffusion. 
As can be seen in Figure 1111 initially the magnetic en- 
ergy increases due to differential rotation and convection 
and only later does exponential decay set in. Depending 
on the field component, the field decay requires nearly 
half of the simulation time to decay back to its original 
amplitude. Therefore, estimating how far the differential 
rotation should have spread if mediated by the magnetic 
field is not straightforward, given the complications of a 
time and space dependent field. In light of the uncer- 
tainty in estimating the timescale of tachcocline spread 
mediated by magnetic field, it is wise to at least investi- 
gate the trend of the magnetic stresses. 

3.5.1. The Lorentz Force in the Tachocline 

The scenario in which a magnetic field aids the spread 
of differential rotation into the radiative interior depends 
on a few factors. First, the field must be sufficiently 
"connected" to the convection zone. Second, the Lorentz 
force must be larger than other forces in the region (at 
least on long time averages). And finally, the configura- 
tion of B^, Bg and Br must be specific. 

As we have seen in Figure [TUl the initial poloidal field 
indeed becomes "connected" to the convection zone and 
is not confined by the meridional circulation. Given this 
clear connection one would expect the differential rota- 
tion of the convection zone to be rapidly imprinted into 
the radiative interior. However, a magnetically-aided 
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tachocline spread depends on the Lorentz force being 
dominant or comparable to the other stresses in the az- 
imuthal momentum equation (at least on long time av- 
erages). This depends not only on the strength of the 
magnetic field in this region but also on the correlation 
of these fluctuating fields. In Figure [9] we showed the 
various forces contributing to the maintenance of the dif- 
ferential rotation. In that figure the Lorentz forces are 
seen as zero because they are substantially smaller than 
the other forces shown. In Figure [Tl]we show the vol- 
ume integrated magnetic energy in the various regions of 
the simulation domain (total: (a), convection zone: (b), 
tachocline: (c) and radiation zone: (d)). We see there 
that the magnetic energy is initially amplified in the con- 
vection zone and tachocline by the action of differential 
rotation and convection. In particular, the energy in the 
toroidal field (initially zero) is amplified to the extent 
that it becomes the dominant field component (dashed 
line). Simply estimating the Lorentz force given these 
amplitudes one would expect it to be comparable to the 
divergence of the Reynolds stress in the tachocline (where 
the velocities are approximately 10 — 10"^ cm/s), but as 
shown in Figure[9l the Lorentz forces are quite small, in- 
dicating that the field fluctuations are correlated in such 
a way that they do not lead to a substantial Lorentz 
force. This can possibly be explained by considering that 
wherever is large it is so because differential rotation 
and/or convection have acted on Br or Bg to generate it, 
so that wherever i?^ is large Br and Bg are not. There- 
fore, while the field clearly "connects" to the tachocline 
and convection zone (see Figure [T0|) , the Lorentz force in 
the tachocline is small and therefore, little angular mo- 
mentum is transported by magnetic terms. This is the 
second possible reason the magnetic field has no effect on 
the propagation of differential rotation into the radiative 
interior. 

The idea that a magnetic field could be connected 
to the convection zone without trans ferring substantial 
angula r momentum was suggested by iGaraud fc Roger j 
()2007f ). Field lines connecting to the convection zone 
are made small scale by convection and turbulence, thus 
imposing an enhanced diffusivity on the magnetic field. 
This enhanced diffusivity could act to disconnect the field 
in the convection zone (in the sense of angular momen- 
tum transfer) from the radiative interior. 

3.5.2. Sense of the Lorentz Force 

Finally, the spread of differential rotation by an inter- 
nal magnetic field depends on the configuration of the 
magnetic field generated at the convective-radiative in- 
terface. Simply looking at the azimuthal component of 
the induction equation (neglecting diffusion): 



nent of the Lorentz force: 



dB^ 
dt 



-rBr-^i^) + 
or r 

d B^ 
or r 



i9Bg d 



r 89 sinO 
sinOug d B, 



i-^n)- (14) 



89 sin9 



one can see that the in the tachocline, where the ra- 
dial gradient of differential rotation is initially dominant 
in inducing toroidal field (first term in Equation 14), a 
toroidal field which switches sign at mid-latitudes is in- 
duced (Figure ITOb ) . Considering the azimuthal compo- 



[(V X B) X B]^ = 



Bg 8{Bfjj sin ( 



Br d{rB^) 
r dr 



(15) 



a toroidal field which switches sign at mid-latitudes, cou- 
pled with a latitudinal field which is positive generates a 
Lorentz force which promotes decelerating angular veloc- 
ity at high latitudes and accelerating angular velocity at 
low latitudes, thereby promoting the spread of the differ- 
ential rotation of the convection zone into the radiative 
interior 

In these simulations we indeed find that the radial gra- 
dient of differential rotation dominates the initial induc- 
tion of the toroidal field, leading to a toroidal field which 
switches sign at mid-latitude (Figure [TUb). Later in 
the simulation, however, the latitudinal component of 
differential rotation becomes dominant. This leads to a 
toroidal field configuration which has one sign toroidal 
field in the northern hemisphere (positive) and one sign 
in the southern hemisphere (negative). Figure llOb . The 
amplitude of the toroidal field peaks at around mid- 
latitudes (see Figure [TOH ) , so that the latitudinal gra- 
dient of toroidal field changes sign at mid-latitudes. Fur- 
thermore, as the field diffuses out of the radiative inte- 
rior and connects with the convection zone, the sign of 
the latitudinal magnetic field changes at high latitudes. 
Whereas initially the latitudinal field was positive ev- 
erywhere in the tachocline, later in the simulation the 
latitudinal field is negative at high latitudes and positive 
at low latitudes. Combining the toroidal field configura- 
tion seen in Figure [TOH . with positive latitudinal field at 
low latitudes and negative latitudinal field at high lati- 
tudes, one can then simply see that this field configura- 
tion gives acceleration at high latitudes and deceleration 
at low latitudes, as seen in Figure [T2l This field con- 
figuration would then act to stop tachocline spread (if 
Lorentz forces were dominant over long timescales). 

Of course, it is not clear that such a field configuration 
would develop in the Sun. But this highly simplified 
simulation indicates that the magnetic field configura- 
tion and associated Lorentz force in a tachocline that is 
bounded on top by convection and fed from below by a 
space and time-dependent magnetic field is quite compli- 
cated and possibly not well described without consider- 
ing instabilities and poloidal field evolution. 

3.6. Weaknesses 

Clearly, this model has several weaknesses. The first is 
that it is essentially two-dimensional and there is some 
possibility that two-dimensional convection conspires to 
limit tachocline spread. It is impossible to test this until 
three-dimensional simulations, using similar parameters 
as those described here are conducted. Such simulations 
are likely not too far off. 

The second major limitation of these simulations is 
that the diffusivities are too large. This prevents us 
from directly testing t he case in which the the rmal hyper- 
diffusion described bv lSpiegel fc Zalinl ()1992j ) dominates 
the viscous diffusion. It also prevents us from getting 
the force balance in the tachocline exactly right. For 

^ The same sense of Lorentz force is found if the induction of 
toroidal field is dominated by the latitudinal gradient of angular 
velocity, because the gradient changes sign with latitude. 
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instance, our meridional circulation appears to be only 
about 45 times to large, leading to Reynolds stresses and 
Coriolis forces in the tachocline which are larger than 
expected. However, this should be more offset again by 
the larger than solar viscous diffusivity. In the end, the 
viscous forces are not dominant in the tachocline in these 
simulations and that is the first order goal. 

Finally, because our timestep is severely limited by 
gravity waves and convection we are not able to run this 
simulation very long. We have only run these simulta- 
tions approximately 20% of a viscous diffusion time or 
40% a magnetic diffusion time. Over this time we would 
have expected more tachocline spread and we are at least 
able to at least the trend of the stresses. The following 
discussion then proceeds with these weaknesses in mind. 



4. DISCUSSION 

These results represent the first attempt at coupling 
the convective and radiative interior of the Sun in order 
to understand the angular momentum evolution of the 
radiative interior. Although troubled by computational 
limitations (as all simulations are), there are still some 



very interesting and important results to be gleaned from 
this model. First, we find that a purely hydrodynami- 
cal model of the solar interior could prevent the spread 
of the tachocline by setting up a counter-rotating cell 
that straddles the base of the convection zone/top of the 
tachocline. The Coriolis force acting on this CRC acts 
to accelerate the flow at high latitude, thus shutting off 
the spread of the tachocline. 

In the magnetic models the meridional circulation ap- 
pears unable to prevent the field from connecting to the 
convection zone; both by diffusion and dredge-up by con- 
vective plumes. However, despite being connected to the 
convection zone very little angular momentum is trans- 
ported between the two regions by magnetic stresses. 
This is because the magnetic stresses in the tachocline are 
small. Despite a rather large toroidal field growing in the 
tachocline due to differential rotation, the radial and lat- 
itudinal field there are quite small and the Maxwell stress 
resulting from the correlations of these field components 
are similarly small. Therefore, while it does not appear 
necessary for a magnetic field to prevent the spread of 
the tachocline, we also find that such a internal poloidal 
field does not promote differential rotation if it connects 
to the convection zone. 
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Friday, July 23, 2010 



Fig. 1. — Model schematic (reproduced from RogersI l|2011l )). Radiation zone occupies the inner 75% of the simulated domain, with 
convection occupying the outer 25%. In the magnetic models a dipolar field is imposed in the radiative interior (field lines shown). 
Gravity-Alfven waves are generated in the radiative interior by impinging plumes. 
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Friday, February 4, 201 1 

Fig. 2. — Time averaged latitudinal flow in non-magnetic (a) and magnetic (b) models. Red represents flow towards the south pole, 
while blue represents flow towards the north pole. The white overlaid line in (a) and (b) represent the base of the convection zone. Flow 
amplitudes are approximately 9 X 10"* cm/s at the top of the convection zone and 5 X 10^ cm/s at the base of the convection zone (colorbar 
units are in cm/s). Bottom four plots show time snapshots from the non-magnetic model with time increasing to the right (c) lO^s, (d) 
lO^s, (e) 9 X lO^s and (f) 1.5 X lO^s. Gravity waves are seen quite clearly in both time snapshots and time averages. 
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Fig. ,3. — Absolute magnitude of the time-averaged radial and latitudinal velocities in non-magnetic (a) and (c) and magnetic (b) and 
(d) models in cm/s. Different line types represent different latitudes; 70° (solid line), 40° (dotted line) and 10° (dashed) lines. We clearly 
see exponential decay of the amplitudes as well as the signs of gravity waves. Amplitudes in magnetic models show slightly faster decay 
and waves have shorter vertical wavelengths. 
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Fig. 4. — Time averaged temperature variation as a function of latitude at Q.004Rq (3 x lO^cm) below the convection zone (solid line), 
O.OliiQ (8 X lO^cm) below the convection zone (dotted line) and O.OSii© (2 x 10® cm) below the convection zone in Kelvin. We clearly see 
warmer poles and cooler equator closer to the convection zone, although there is substantial variation at low latitudes. Further from the 
convection zone this profile reverses, with slightly cooler poles and warmer equator. 
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Friday, February 4, 201 1 

Fig. 5. — Time averaged azimuthal flow in non-magnetic (a) and magnetic (b) models. Red represents positive values with respect to 
the mean Q (441nHz), while blue represents negative values with respect to the mean. The overlaid white lines represent the base of the 
convection zone and the base of the tachocline after it has stopped spreading, or O.14R0 below the convection zone. The bottom four plots 
show time snapshots from the non-magnetic model with time increasing to the right (c) lO^s, (d) lO^s, (e) 9 X lO'^s and (f) 1.5 X 10*s. 
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Radius R/Rq 

Fig. 6. — Absolute value of the radial gradient of the angular velocity as a function of radius beneath the convection zone for the non- 
magnetic model (solid line) and the magnetic model (dotted line), in units of \/cm/s. The low latitude values are averaged from 60°-120° 
and the high latitude values are averaged from 0°-30° and 150°-180°. The half-maximum value is defined as half of the amplitude at the 
secondary pealc (denoted by the dashed line) and marked with the solid line. Defining the tachocline depth as distance between the base 
of the convection zone and this position, the tachocline penetration depth is measured as 0.021iiQ for high latitudes and O.O22il0 for low 
latitudes in the non-magnetic model. 
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Fig. 7. — Tachoclinc thickness as a function of time for the non-magnetic model (soUd Unc) and the strong magnetic case (dashed hne) 
in cm. This thickness is measured as described in the text and is averaged over the polar regions (0°-30°, a) and the equatorial regions 
(60°-90°, c). The time derivative of the tachoclinc penetration depth is shown in panels (b) and (d) in cm/s. Polar regions show that the 
tachoclinc thickness has stopped increasing. Lower latitudes still show a slow increase. The insets in panels (b) and (d) show the time 
derivative of the tachocline penetration depth over just the last half of the simulation. 
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Fig. 8. — Anisotropy of the Reynolds stress, measured as the ratio of the first two terms in Equation 12 to the last two terms in Equation 
12. Different line types denote different radii within the ta<;hocline, O.71il0 (solid), O.TORq (dashed) and 0.69Rq (dotted). 
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Fk;. 9. Forces within the tachocHne in the non-magnetic model (left-hand panels) and in the magnetic model (right hand panels). 
Three difFerent radii within the tachocliiie are shown. Solid line represents the total force, i.e. the llis of Equation 13 (all forces are shown 
in units of grn / {crns)^) . Dashed line represents the Coriolis force (the sum of the first and fourtfi terms on the rhs of Equation 13). Dotted 
line represents the Reynolds stress (second term on rhs of Equation 13) and triple-dot dashed line represents the diffusive terms (the sum 
of the fifth and sixth terms on rhs of Equation 13). In the magnetic model the Maxwell stresses (third term on rhs of Equation 13) axe 
also shown (although the appear as zero). 
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Fig. 10. — Evolution of toroidal (top panels) and poloidal (bottom panels) field. The dark line represents the base of the convection zone. 
Toroidal field evolution: Initially, there is no toroidal field (a). Differential rotation acts on the poloidal field to produce a toroidal field 
which initially changes sign at mid-latitudes (b). This configuration then becomes unstable and leaves behind a single dominant sign in 
each hemisphere (c). The time averaged toroidal field is shown in (d). Poloidal field evolution. After a short time the poloidal field diffuses 
into the convection zone (f), aided by convective overshoot and dredge-up, thus connecting the radiative interior to the convection zone. 



Solar Rotation 



19 




Time (x 1 O^^s) 

Fic;. 11. Magnetic Energy as a function of time for the three different field components; radial field (solid line), latitudinal field (dotted 
line) and toroidal field (dashed line). The total integrated energy for the entire domain are shown in (a), the convection zone (b), the 
tachocline (c) and the radiative interior (d). One can clearly see an initial increase in all components of the magnetic field energy in the 
convection zone and tachocline followed by exponential decay. The toroidal field is amplified substantially and shows oscillatory behavior 
where the toroidal energy grows, decays and repeats. In the raditiative interior, the field predominantly decays, except the toroidal field. 
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Fig. 12. — Lorentz forces within the tachochne. SoUd Hne represents the total Lorentz force (rhs of Equation 15), dashed hne represents 
the second term on the rhs of Equation 15 and the triple-dot dashed line represents the first term on the rhs of Equation 15. The dotted 
lines denote zero. The left two panels show high latitudes, while the right two panels show only the lower latitudes. The top two panels 
are at the top of the tachocline, O.71il0, while the bottom two panels show the components of the Lorentz force at the bottom of the 
tachocline at O.69R0. 



